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abstract 

The stability has been investigated of the unsteady flow past an infinite flat plate 
when it is moved impulsively from rest, in its own plane. For small times the in- 
stantaneous stability of the flow depends on the linearised equations of motion which 
reduce in this problem to the Orr-Sonmierfeld equation. It is known that the flow for 
rtam values of Reynolds number, frequency and wavenumber is unstable to Tollmien- 
Schlichting waves, as in the case of the Blasius boundary layer flow past a flat plate 
w.th increase in time, the unstable waves only undergo growth for a finite time inter- 
val, and this growth rate is itself a function of time. The influence of finite amplitude 
effects is studied by solving the full Navier-Stokes equations. It is found that the 
stability characteristics are markedly changed both by the consideration of the time 
evolution of the flow, and b y the introduction of finite amplitude effects. 
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§1 Introduction 

In the flow over a flat plate at zero pressure gradient the stability of the boundary 
layer is not only dependent on its instantaneous thickness, but also on its rate of 
growth. This problem, based on the Blasius mean flow has been studied extensively. 
However, experiment shows that the stability characteristics based on the local growth 
of the boundary layer provide a good approximate model. This was confirmed in the 
triple-deck structure calculations of Smith (1979a) on the influence of the rate of 
growth of near neutral unstable Tollmien-Schlichting (TS) modes at high Reynolds 
numbers. The weakly nonlinear evolution of TS modes within a Blasius layer was 
studied by Smith (1979b) and Hall & Smith (1984). In the latter of these studies, it 
was demonstrated how a linear mode could evolve smoothly into a finite amplitude 
perturbation. In Otto (1994), this was shown to be feasible within the unsteady 
problem. Smith also found that it was not possible to determine the effects of boundary 
layer spatial growth on the stability of the modes since different conclusions were 
found depending on which flow quantity was used to express the finite amplitude of 
a mode. Gaster (1974) and Eagles & Weissman (1975) have also contributed to our 
understanding of the Blasius problem. In this flow, Tollmien-Schlichting waves grow 
spatially in the flow direction and the resulting equations in this coordinate are elliptic, 
and are computationally expensive to solve. The approach of Bertolotti (1991), using 
the Parabolised Stability Equations (PSE), in which the elliptic equations are made 
parabolic on the assumption of slowly evolving wave characteristics and the changes in 
structure may be locally assumed linear, has been one of the more successful methods 
to reduce the computational cost. This approach allows far larger steps to be taken in 
the evolutionary coordinate. 

In the present work, the related temporal problem is considered. Otto (1994) found 
reasonable agreement between triple-deck and Orr-Sommerfeld results in regimes of 
common validity, for the unsteady flow problem of the flat plate moved impulsively 
from rest. Here, we discuss both the linear and nonlinear problems, and since our 
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problem is already parabolic in tbe evolutionary coordinate, we use the PSE method 
ology to study the effects of finite amplitude and to reduce the number of time steps 
per period required for accuracy, which it does very effectively. We find substantial 
differences between the results of linear theory based on the Orr-Sommerfeld equations 
and our results for finite amplitude waves. 

In Section 2, we derive equations for the linear and nonlinear evolution of Tollmien- 
Schlichting waves with a brief summary of our period fitting technique. Section 3 
contains details of the numerical methods used to solve the equations. In Section 4, 
we discuss the results and finally in Section 5, we draw some conclusions. 
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§2 Formulation 


The physical problem consists of an infinite flat plate in two dimensions which is 

at rest within a fluid of density «, and pressure p„. At time t = 0, the plate is moved 

impulsively from rest with the speed £/., in the direction parallel to its alignment (see 
Figure 2.1). 


Po 

Po 


//// / / / / / / /// 7 - ~ Q> > 

Figure 2. 1 Velocity Profile Due to Impulsive Motion of Flat Plate 


The governing equations for this problem are the equation of continuity and the 
Navier- Stokes equations 


du do 
dx + dy ^ 
du du du 

dt dx dy 


dp 1 

o ^ u 

dx Re 


( 2 . 1 ) 


dv dv dv dp 1 

where each variable is dimensionless, with sealing parameters V „ p„, and L repre- 
senting the ambient density, plate velocity, viscosity and a unit of length in the flow 

direction, respectively. The Reynolds number per unit length of plate, Re is defined 
as 

LUopo 


Re = 


Vo 


( 2 . 2 ) 


and V 2 is the two-dimensional Laplacian operator 


V 2 = 


d 2 d 2 
4 ~ 


dx 2 "h dy 2 ’ (2-3) 

We note that the pressure has been non-dimensionalised by Po U 2 and time by 
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I„ the case of an infinite plate impulsively set in motion, the mean solution, «. 
the undisturbed flow, is given in dimensionless terms by the reduced set of equations 

du _ 1 
dt _ Re dy 2 
v = 0 

P o 


(2.4) 


P = 


with initial and boundary conditions 


u 


( 0,0 = 


PoU% 


0 t = 0 

1 t > 0 


(2.5) 


The mean streamwise 


u(Y,t)~* 0 T-^oo 

se velocity component is therefore given exactly by 


, xiyRt \ 

( y ’ *) = 1 - erf ^ _ 2VT 


( 2 . 6 ) 


Given the form of (2.6), it is appropriate to choose physical variables X and Y that 
absorb the \fRe dependence, 


X = xsfR~e 

Y = yVRe 

With this transformation, the governing set of equations become 

du dv 


(2.7) 


+ — = 0 
dX dY 


1 du 


du 


du 


dp 


■jRidt dX 

+ U 


— + dX + JTe 

1 dv dv dv dp 1 




(2.8) 


y/Re dt 

and the mean flow is now given by 


d X +V dY-~eY + Vlfe 


ii (y, t ) = l-erf(T f ). < 2 - 9 ’ 

The lack of a physical boundary condition for pressure provides a computational 
difficulty for equations (2.8) cast in terms of the primitive variables. We avoid this 
problem by reducing this set of equations to a single equation for the normal velocity 
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component. In this section, we will discuss the analysis for both linear and nonlinear 
disturbances. In addition, we show how to generate an appropriate initial condition 
to begin the stability calculation, and introduce a period-fitting technique intended to 
reduce the total computational effort. 


2.1 Linear Disturbance Equations 

We let v and p represent disturbances from the mean state: 


u(X i Y,t) = u(Y,t) + ii(X,Y,t), 
v(X,YJ)=v(X,YJ), 
p(X,Y,t) = J^+p(X,Y,t). 


( 2 . 10 ) 


After substituting equations (2.10) into equations (2.8) and linearizing the resulting 
expressions, the perturbation pressure, p. and the streamwise velocity component, u, 

may be eliminated, giving the following single equation for v, the transverse compo- 
nent: 


/ 1 d d 
' \/Re dt + U dX 


, 3 2 N d 2 d 2 

s/Re. ' dX 2 + dY 2 ') ‘ dx 2 + dY 2 ' 


d 2 u d 
dY^dX 


( 2 . 11 ) 


The physical problem is spatially homogeneous in the X-direction, and it is legitimate 
to assume a periodic X- variation, giving 


v(X,Y,t) = v(Y,t)e iaX 


(2.12) 


which reduces (2.11) to 


'( 


i d 
VRedt 


au 


\( 92 2 \ 9 ’ 2 U 

)(^p2 ~ a ) + a 


+ 


dX 2 iVR~e K dY* 


d 4 


- 2a 2 


d 2 

dY 2 


+ 0 4 ) 


v = 0. (2.13) 


When a periodic time dependence, e - <w< , is assumed, and noting the mean velocity 
is a function of Y and t only, equation (2.13) reduces to the familiar Orr-Sommerfeld 
equation with similar results. In our case, we consider the solution to equation (2.13) 
at any fixed time and then the results show the local stability characteristics of the 



flow for the linear problem. For this time-dependent problem, the boundary conditions 


are 


5(0, t) = 0, 


(2.14) 


and since ^ = 0 at Y — 0, 


dv 

ay 


-(0,<) = o. 


(2.15) 


In the farfield, as Y oo, the solution depends on e- 12 , giving for large Y, 

Y — > oc . (2.16) 


+ a $(Y,t) = 0 


av 

a 2 »(y,t) 

a 2 v 


cv 2 6 


(Y,i) = 0 


2.2 Nonlinear Disturbance Equations 

Substituting (2.10) into (2.8), retaining the nonlinear terms, and eliminating the 
perturbation pressure, we find the following for v. 


[/ i d _ d _ 1 /_g . d 2 a L + 

[(yUa? + “av ar2 ' ajf2 

, a 8 u o 2 d 2 \ . ae / a 2 

‘ ( s ax + 5 aid fax* + aF* r + ar W 


a 2 a*s _a_i _ 
M 72 _ ar 2 J v 

d 2 v . 2 ^fjg_ . 
ar 2 j u+ axUx 2 


a 2 


ar 2 


u . 


(2.17) 


The g* could be eliminated using the continuity equation. However, since the ** 
term cannot, we leave both terms in for ease of presentation. 

We now express the disturbance quantities as a finite sum of modes with respect 


to the wavenumber alpha: 


M 

v(X,Y,t) = ^ ij(Y,t)e iiaX , 

M 

u(X,Y,t) = 5Z U){Yj)e vaX , ( 2 - 18 ) 

j=-M 

M 

p(X,Y,t) = 52 p,(Y,t)e'’" x . 

j=-M 


The continuity equation allows for elimination of all the u/s, except u„, in favour of 
the Oi’s and requires that S„(F,() = 0. The resulting equation after making these 
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substitutions in (2.17) and neglecting all products which combine to give a mode 
beyond M is given by 


/ i d 

\y/Rtdt 

where 


— jau + 




d 2 


■2 2> 
- J Oi 


))(^2 - J 2a ' 2 ) +J 


a 


d 2 u 

dY 2 


Vi = 


= Nv (2.19) 


M+j 


N v = i 


k = ~ M 
k^0,j 


•i ifdv. 
jka 


dY 


* „ k dvi~k\ 

7Vj-k — : — rVk —^rr~ ) + 


j d 2 v k dv } - k j d 3 


Vk 


3 - k dY J j - k dY 2 dY k dY 3 


i ' / ~ & i>j ■22-'' 0“Uq ^ \ 

+ ja \ Uo dY^ ~ J a UoVj ~ dY^ Vj ) 


d 2 i 


v i-k 


( 2 . 20 ) 


for j = — M, -M + 1, -M + 2, -1, and 

M 

I i km 2 ( 

dY 


k = -M+j 

k*0,j 


Nfi = i y j kcx 2 ( -Tr-Vi I - l JdH k 

Var J k j-k Vk dY ) + j-kdY 2 dY ~kdY 3 


Vk „ 
J v J-k 


d 2 dj 


d 2 u 0 


+ gyi - J 2 “ 2 « oij ~ ^0,) 


( 2 . 21 ) 


for j = 1,2,3, ..., M. 

The boundary conditions for the t)/s are analagous to those for the linear case: 

«i(o,0 = o, |lf-(o,<) = o, 


Svj(Y,t) 

dY 

d 2 vj(Y,t) 

d*Y 


+ javj(Y,t) = 0 

- j 2 a 2 Vj(Y,t) = 0 


OO . 


( 2 . 22 ) 


The vaiiable v oCY,t ) has no associated wavy part, and is effectively a correction to 
the mean solution due to the presence of the finite amplitude waves. An equation for 
it may be found by looking at the non- wavy component of the ^-momentum equation, 
after substituting in (2.18): 


du 0 d 2 u 0 
~dt ~ dY 2 


M 

— y/We. 


Jfc = - M 
/t^O 


i d 2 v k „ 
ka dY 2 V ~ k 


(2.23) 
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The boundary conditions for u 0 consist of a no-slip condition at the wall and an 
asymptotic condition far away from the wall. 


u o((M) = 0, 


dupjYj) _ 
dY 


Y 


oo . 


(2.24) 


The system of equations (2.19) through (2.24) define the nonlinear problem. 


2.3 Initial Condition 

In order to solve either the linear or nonlinear problem for some time t = t 0 > 0, 
we need to define an initial perturbation which is at least an approximate solution 
of the linearised equations. To this end, we can freeze the flow at a given time and 
determine the local stability characteristics of the flow from an Orr-Sommerfeld (OS) 
analysis. In particular, for given values of « and >/He, we can determine a function 
v(Y) and a frequency w, such that an initial disturbance of the form u(K)e lc * x will 

behave according to for a11 time for the fr ° Zen fl ° W ‘ 

However, since the OS equation does not allow for the evolution of the undisturbed 

flow with time, using an initial condition constructed in this way generates spurious 
oscillations in the solution. We overcome this problem by finding an initial condition, 
which is a slight perturbation of the OS eigenfunction v(Y), that is sensitive to the 
changing mean flow. Following the approach of Bertolotti, we let 

v{Y,i) = v'(Y,t)e~'fo d(r,)dri (2.25) 


where it is assumed that the time dependence of v 1 is slowly varying. Substituting this 
expression into the linear stability equation (2.13) gives 


(m_ 


- au(Y,t)j ( 


d 2 

ri 


— a I v 


\Y,t) + a§£v'(Y,t) 


, _J —(* 


dY 2 / 

-2 cv 2 ^- + cv 4 )u(^<) = -^(wi 


(2.26) 
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We now make a I st order Taylor Series approximation about t 0 : 


t ~ ^0 + T 


0 ^ 0(to) + r (tfo) = + t$i , 

v'(Y,t) « v'(Y,to) + r^(Y,to) = v' 0 (Y) + rv[(Y ), 
™ %-(Y,t 0 ) = v{(Y) , 
u(F,t)«t2(F,* 0 ) + rjf(F,t 0 ), 


d 2 u(Y,t) 

dY 2 


d 2 u 
8Y 2 1 


k \ ^0 ) + T 7 


a 3 u 


1*4 2 


( V, tn ) 


(2.27) 


7) t / y • 2 f 

We note that the approximation for ■ does not include a ^-(F,i 0 ) term; the 

slow variation of v f with t allows us to ignore the second derivative. Collecting like 
powers of r, namely r° and r 1 , we have the following equations, respectively: 

' ' (2 98) 

where L is the Orr-Sommerfeld operator 


L = 



— ou(F, <o 



+ a|j 4 (r,*o) + 



d 4 

dY 4 


2 a 


2 


a 2 

dY 2 


+ o 4 ). (2.29) 


The boundary conditions for and v\ are given by (2.14), (2.15), and (2.16). 

With the additional constraint that Vq and v\ are orthogonal, 

/*oo 

/ (v£)*i/JdF = 0, (2.30) 

Jo 

equations (2.28) provide an eigenvalue problem which, when solved, gives an initial 
profile v' 0 (Y) and frequency 0 0 which, to first order, take into account the evolution of 
the undisturbed flow. 

To further improve the quality of the initial condition, a higher order version of 
the above can be performed. In this case, we use a second order Taylor Series approx- 
imation, and include the second derivative of v 1 with respect to f, but ignore the third 
derivative. In addition to 6 0 , 0] , v ' 0 and v \ , the resulting eigenvalue problem contains 
unknowns 62 and v 2 which represent the second-order Taylor Series terms. Extending 
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the above analysis to 2 nd order gives a linear equation for v' 2 and an orthogonality 
condition, analagous to (2.30). 


2.4 Period Fitting 

As mentioned in the previous section, we can separate v(Y,t), for the linear case, 
into a time periodic part, and a function v'(Y,t) that varies slowly with time. In this 
case, we will write v as 

v(Y\t) = v'(Y,t)e~ luJt (2.31) 


where u is taken to be a real quantity; the growth or decay is represented by the 
t - variation of v' . To adequately resolve the oscillatory part requires many time steps 
per period of oscillation. However, the temporal resolution requirements of v' are 
far less restrictive. Consequently, it is more efficient to remove the periodic part of 
the solution from the equations and compute only the slowly varying function v' . 
Substituting (2.31) into the linear disturbance equation (2.13), we find the following 
equation for v' 


v^V ay2 



dv'(Y,t) 

at 


+ (-fa - cu(Y, ()) (<#= - n 2 )»'(r.<) + 


+ 


i ( a 4 

i^FU\dY< 


-2a 2 £?+a 4 )v'(Y,t) = 0. 


(2.32) 


Since the extent of the flow increases as y/t, the value of u changes with time. 
Therefore, to implement (2.32), we must be able to estimate u> at any given time. We 


define the energy E(Y,t) in the following way: 


roo r°° 

E(Y,t) = f v\Y,t)4Y= v'\Y,t) e - 2i -‘AY . (2,33) 

Jo Jo 

Ignoring the slow changes of v f and u> with time, we approximate u; by — 2 iE ~dt * 
However, we wish now to work only with v . Initially, we use a value of to in (2.32) 
given by the OS analysis or by the method described in Section 2.3 above. At each 
time step, we can correct this value by recognizing that the computed value of v 
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will contain the slowly changing envelope plus an oscillatory part with a very small 
frequency that represents the change in u i.e. 


f°° 

E'(Y,t)= / v' 2 (Y,t)dY , 

Jo 


Alo 


i be* 


2iE ' dt ’ 
^new = ^old H“ ■ 


(2-34) 


For the nonlinear case, the procedure is very similar. Given that we have 


ot } = jo 


(2.35) 


we make the assumption 


LOj = ju). 


(2.36) 


That is, we compute uq , associated with the fundamental mode, by the method shown 
in equation (2.34). Then all uq for higher order modes are given by (2.36). We note 
that any errors incurred by this approximation are at least partially corrected by the 
the time dependence of t>'. 

The nonlinear disturbance equations for the u's now become 


l 

\/ Re 


(w* -)V)®5 +}(^ -cm(Y,t))(g, - ]V)v'(Y,t) + j 

= N.,' 


a 2 s(v,«) 

ay 2 


v'(Y,t) 


i ( a 4 

;\/k7 V 5V " 1 


(2.37) 

where Nt/ has the same form as given in (2.20) except that v is replaced everywhere 
by v' . 
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§3 Numerical Methods 

In this section, we will provide some of the pertinent details of the numerical 
solution of both the linear and nonlinear problems. We begin with a brief discussion of 
the grid and spectral discretization, and follow with descriptions of the initial condition 
generation, linear simulations and nonlinear simulations. 

3.1 Discretization 

The mean solution is given by equation (2.9) which is repeated here: 

u(Y,t) = l-etf(~i). (3.1) 

Although, u is technically finite for any finite Y and t, it decays quickly away from the 
plate, and from a computational standpoint, it is appropriate to define the extent of 
the undisturbed flow to be the Y value at which the undisturbed velocity is 1% of the 
plate speed. Since typically the OS eigensolutions decay slowly relative to the decay of 
the undisturbed flow solution, it is necessary to choose a physical domain that extends 
well beyond this limit. Defining the edge of the undisturbed flow at t = t 0 to be 

= i - = o-oi. ( 3 - 2 > 

we have F ed = 1.821(2 v / to). We choose an upper bound for Y to be 

Y m ax = 40\/to, (3.3) 

i.e. our domain extends to about 10 times the initial undisturbed flow extent. 

We employ a Chebyshev collocation discretization using N — SO polynomials. Al- 
though the disturbances do not require any local resolution because, once introduced, 
they extend throughout the domain, the mean solution decays rapidly away from the 
wall. Consequently, we use a mapping which maps 0 to l^max onto ^ required 

by the spectral discretization, but does so in such a way that points are concentrated 
near the lower boundary. The computational coordinate ?/ is defined as follows: 


where b and c are chosen so that F(l) = F max and F(0) = F mid The value of F mi d is 
the physical coordinate that bisects the grid. We take it to be 15y/to, or about 4 times 
the initial undisturbed flow extent. This allows for the undisturbed flow to grow in 
extent during the simulation, and keeps the stretching from being too severe. With 
these stipulations we have 


c = 


bn ax bn id 

^max 2 Vn-jid 


b — 1 + 2 — 

'max 


(3.5) 


0 


H h 


Y 


mid 


1 — H+ 


Y 


Y 

max 


Figure 3. 1 Representation of Stretched Spectral Grid 


3.2 Initial Conditions 

For both the linear and nonlinear simulations, the disturbances are introduced at 
dimensionless time to = 0.25. We generate three sets of initial conditions at this t 0 - 
The first corresponds to the OS eigenfunction, and the remaining two are the l s< and 
2 nd order corrections, respectively, of §2.3. 

The continuous Orr-Sommerfeld eigenvalue problem is defined by 


L v'o(Y) =£ 0 


(3.6) 


and the boundary conditions 


f«0) = 0, §£(0) = 0, 


3 2 t4(V) 

dY* 


0 


(3.7) 


oo 


a 2 v' 0 (Y) = 0 

where L is the Orr-Sommerfeld operator given by equation (2.29). The corresponding 
discrete problem replaces ^7 with the discrete operator (^ 7 ),, where given a discrete 
function <j> defined on the nodes / = 0, 1, 2, ..., iV, 

(d y) 1 = idy) (3.8) 

k = 0 
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where, once again, ?/ € [-1,1] is the computational coordinate. The d/,*’s are the 
coefficients of the Chebyshev 1 st derivative matrix. The system is then defined by 
applying the discrete OS equation at the interior nodes 2, 3,4, ..., N — 2, and the discrete 
wall and asymptotic boundary conditions at t] = — 1 (Y = 0) and r/ = 1 (Y = V max ) 
respectively. 

For a given a and y/Re, the solution of this eigenvalue problem provides a frequency 
u> and an initial profile v°,Z = 0, 1, 2, ..., N, which may be used to start the linear- 
simulation. 

In the nonlinear case, this profile may be used as the starting value for the funda- 
mental mode. In these simulations, we set the higher order modes initially to 0, and 
the —1 mode to the complex conjugate of the fundamental. 

The discrete analogues to the 1 st and 2 nd order correction problems of §2.3 follow 
directly from above. The resulting linear systems are also eigenvalue problems, the 
solutions of which provide a frequency and initial profile. 


3.3 Linear Simulations 


The continuous linear problem is defined by equation (2.13) with boundary condi- 
tions (2.14), (2.15), and (2.16) and initial condition provided as in the the above sec- 
tion. We employ the discrete form of these equations along with the Crank-Nicolson 
time differencing, given by 


-n-f 1 
"i - 

At 


1 ( di i ' i 

2 V dt ' dt /’ 


n + 1 


di) V 


(3.9) 


which results in the following linear set of equations for £" +1 : 

(7E7 “ 2 CVU " +1 ^) ((m 7 ) “ a2 ) + 2 ft (fv^)/ ^ 


+ 2 t yfe((D D vO 


2< * 2 (w) + « 4 )ju" +1 
= [(yfe + 2 ftU "^) ((w) “ ft2 ) “ 2 a (w^)l ^ 


(3.10) 


At 


2i \/Re 




Vl 


l = 2 , 3 , 4 ,..., N - 2 , 


■* riH- 1 


and 


0, 


V DV' )o 


(* + n)«r= 0 . ((^) 2 -« 2 )«« +, = 0 . 
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The value of N for all simulations is chosen to be N = 80. 


To perform the simulation with period fitting, we implement the discrete form 
of (2.32) along with the Crank-Nicolson time differencing (3.9) and the appropriate 
boundary conditions. This linear set of equations for (v')" +1 , / = 0, 1,2, ..., N is given 
by 


(yk + “ a “/ n+1 )) ((i &) 2 - a 2 ) + |a(|^-)” +, Af 


'\n + 1 


/ / \ 7i 

( V )l 


and 


_2or2 (rfp) + ° 4 ) 

(yk ” ^(yk - a “")) ((^) 2 - « 2 ) - ^o'(|^4)"At 

- 2n7k((rfy) “ 2Q;2 (m 7 ) 2 + « 4 )J( u ')r l = 2, 3, 4,..., TV - 2, 

(Oo +1 =0, ^K)o +1 =0, 

+ o)(u')^ +1 = 0, ((jjy) — a 2 j(t/)) 


(3.11) 


)* + 1 =o. 


3.4 Nonlinear Simulations 

The continuous equations for the nonlinear case are given by equations (2.19) 
through (2.24). The time differencing will once again be of the form (3.9). The result- 
ing discrete set of equations is nonlinear, requiring iterative solution. Consequently, 
we will verify results by also performing a simulation which uses Crank-Nicolson for 
the linear left-hand side of equation (2.19) and explicit 2 nd order Adams-Bashforth for 
the nonlinear terms on the right-hand side. 

For the explicit treatment of the nonlinear terms, we have the following discrete 
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analogue to equations (2.19) through (2.21): 

- Ijasr-At) ((*)’ - jW) + ya(&);+'At 

+^((w) 4 - 2 > 2 “ 2 (^) 2 +i 4 « 4 )]*;r 

= (75? + jj““" A( ) ((ot) 2 _ j2 “ 2 ) “ 5J“(li*)l A( 
-i^((w) , - 2 J 2 “ 2 (^) 2 +^ 4 ). 

for / = 2, 3, ..., N - 2, where 


(3.12) 




M + ; r 


NO", 


Dt5£, />n 

],1 ~ \J -- \ D y ' v )-k,l - j _ k u k,l D y 


= * X! j fco;2 ( 


* r D ”" 

v kA 


k = — M L 
fc^O,; 


, j pa *g.i P*i-*.* 

3F ) + j — k D F 2 DF 


jD^I-n 

fc DF 3 


/ D 2 0", , 2 2 „ „ B D 2 up , w \ 

+ jayu 0 jyy2 •? 01 UqV },‘ Y)Y2 V hl)’ 


(3.13) 


for J = -M , -M + 1, -M + 2, ..., -1, and 

M 


NO" 




fcan-Af+i 

**0,j 


DF Uj - u j-k 


DO?* A j 

D y ; + j_ 

• A -2 2- * n ° 2{i 0-„ 

+ }(*{ 


j 

fc DF 2 DF 


n3. A , n 

J U V k,l ~ n 

k DF 3 j_M 


D 2 0?, -2 2- /. n D^ 0 .„\ 

“°-Dy§" ~ 3 a UoV i’ 1 ~ DF^')’ 


(3.14) 


f or j — 1,2. 3. .... A/. In these expressions, the superscript represents the time level, 
the first subscript indicates the mode, and the second subscript refers to the grid. 
Similarly, the discrete values of the mean flow correction are governed by 


(l - ^(w) 2 )“oT _ (} + ^ ( DV ) ) U 0,l 

i\/7teAt V /3-n P 2 </ 1 .V»-l 

25 — Z, \J u (-fc,0 dy 2 fc u (-M) DY 2 J 


(3.15) 


k = -M 
k^O 
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The boundary conditions for each mode j E [- M,-M + 1 are given by the 

discrete forms of (2.22) for v" +1 


(£2!S 

V DY / 


-n+1 

J j, 0 


( D£+> 


N + S' 


= 0, (■ 


DY 2 


DY 
) /V — a 


)o= 0 ’ 


2-n+l 


= 0, 


(3.16) 


and (2.24) for the mean flow correction: 


u 


71+ 1 

0,0 


o - 0 . 

V DY ) N U 


(3.17) 


The above equations constitute a linear set of equations for the unknowns i)"* 1 , j = 
— M,—M + 1 , / = 0, 1,2,..., TV, where, again, TV is taken to be 80. We use this 

simulation only as a check. The period-fitting technique discussed in §2.4 allows for 
a drastic increase in time step. However, stability considerations will not allow full 
advantage of this increase to be taken when the nonlinear terms are treated explicitly. 
The majority of the results presented in the following Results section were obtained 
using period fitting and Crank-Nicolson differencing for the nonlinear terms, resulting 
in a fully nonlinear set of algebraic equations to be solved at each time step. In this 
case, (3.12) is replaced by the discrete analogue of (2.37) 



2iVfc((l>v) ~ 2 J 2a2 (wy) +/<* 4 ) (»')”,/ 


+ f (N(u')r + N(u')r) , 


(3.18) 


and (3.15) is replaced by 


(i -¥(&)>?,!' = (i+¥(&) 2 ) 


u n 


0,1 


M 

i^At V- J_ ( *n D ^h 
2a 2—4 k \ V (-k,l) DV 2 
k = ~ M 
kytO 


+ v 1l+1 - 


k ,1 


DY 2 


-)■ 


(3.19) 
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To solve the nonlinear set, at each time step we begin by making a guess for the 
solution which is a cubic extrapolation using the previous 4 time steps. Using this 
guess, we compute the nonlinear terms on the right-hand sides of equations (3.18) and 
(3.19), and then solve the resulting linear systems for new values of u" , ' which serve 
as the new guesses in the iteration. The process is repeated until convergence. 

To utilise period-fitting, we employ the discrete analogue of equation (2.37). The 
algorithm for solving the resulting nonlinear algebraic set of equations at each time 
step is then identical to that described above. 
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§4 Results 


In this section, we present results for both the linear and nonlinear simulations. 
We begin, however, by presenting the initial conditions obtained by Orr-Sommerfeld 
analysis, and by the 1 st and 2 nd order corrections. 

4*1 Initial Conditions 

The discrete version of the Orr-Sommerfeld equation (3.6) and corresponding bound- 
ary conditions (3.7) is solved by fixing values of a and y/Re and solving the eigenvalue 
problem for the frequency uj and eigenfunction v' 0 . The value of u; is generally complex 
with a positive imaginary part representing growth with time and a negative imagi- 
nary part indicating temporal decay. For particular pairs of a and y/Re a purely real 
lj results. Figure (4.1) shows the neutral Orr-Sommerfeld curve for time to — 0.25. 
The results are identical to those given in Otto (1994). We present this figure for later 
reference. 



Re 6 


Figure 4.1 Neutral Orr-Sommerfeld curve showing the pairs 
of Reynolds number and streamwise wave num- 
ber which give purely real a> 
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The reference values of a and \fRe we have chosen for our simulations are 


a = 0.4 
yfRe = 3000 


(4.1) 


which defines a point just to the left of the neutral curve in figure (4.1). For these 
values, the corresponding u) resulting from the OS solution and the I s * and 2 order 
corrections are given by 

Orr- Sommerfeld: 0.27589729 — 0.00004782? 

l a< Order Correction 0.27585418 + 0.00028591? 

2 nd Order Correction 0.27584836 + 0.00028872? 

and the resulting eigenfunctions, which may be used as initial conditions to the linear 
and nonlinear simulations, are presented in figures (4.2a) and (4.2b). Figure (4.2a) 
shows the real part of the initial profile, as generated by the OS analysis, which we 
call v$ s . In figure (4.2b), two curves are given. The dashed curve shows the deviation 
of the real part of the initial profile generated by the 1 st order correction, v ( ] : from 
that generated by the OS analysis, u° s . The dotted curve in this figure shows the 
deviation of the real part of the initial profile generated by the 2 order correction, 
Uq, from that given by the I st order correction. To plot both of these on the same 
graph, we scale the two curves appropriately, so the actual quantities seen on this 
figure are |10 3 I and 1 1 0 4 j vs. Y, where the denominators represent the 

peaks of these curves. This figure shows that the 1 J< order correction represents a 
0.1% deviation from the OS curve, and the 2 nd order correction deviates from the I s ' 
order curve by about 0.005%. It will become evident in the linear simulations that the 
subtle differences between the initial conditions generated by each of these methods 
are important. It should also be noted that the OS equation predicts a stable mode 
while the corrected value represents an unstable one, for this pair of a and \/Re. 
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1.0 


0.8 
Deviation 


Figure 4.2 Initial Profiles for v. (a) Generated by OS analysis. 

(b) Scaled deviation between 1st order correction and 
OS (dashed), and between 2nd order correction and 1st 
order correction (dotted). 




4.2 Linear Simulations 


The linear simulations without period fitting are solutions to equations (3.10) at 
each time step. As results, we present the growth rate, Im(u;), as a function of time, 
where to is computed at each time step by equation (2.34). We note that a positive 
value of Im(u;) indicates instability. We also point out that if u is fixed at its initial 
profile during the simulation, and not allowed to grow, this is equivalent to the parallel 
flow assumption and the value of to should be constant in time and equal to the OS 
value. This provides an initial check of the simulation code, and is what we observed. 

For the simulations without period fitting, the dispersive error in the Crank- 
Nicolson time differencing scheme makes the results very sensitive to the number of 
time steps used per period of oscillation. We begin by presenting results for a simu- 
lation which uses 200 time steps per initial period of oscillation; i.e., At is chosen so 
that based upon the value of to at t = to, one period of oscillation is divided into 200 
grid points: 

^ = 200Real(u/ 0 ) * (^-2) 

For this case, we give numerical solutions using each of the three initial conditions 
detailed in §2.3. Figure (4.3a) gives the results for the initial condition generated 
by the OS analysis (solid line), and the 1 st order correction. It is evident that the 
nonphysical oscillations in the OS case are a result of the poor starting condition. The 
l 5 * order correction, which is sensitive to the evolution of the mean flow, eliminates 
most of the oscillations. The simulation using the initial condition generated by the 
2 nd order correction is free of oscillations on the graphical scale. This result is given 
in figure (4.3b). 

Figure (4.4) gives the solutions corresponding to 100 (solid) and 50 (dashed) time 
steps per initial period of oscillation, respectively. The dotted curve on this figure is 
the 200 time step per period result given in figure (4.3b). We see that using 100 time 
steps per period is still ample, but that the quality of the results when 50 are used is 
markedly diminished due to dispersive error in the time-differencing scheme. 


22 



In figure (4.5), we give the results for the linear simulations with period fitting, 
using a time step equivalent to 10 (solid) and 2 (dashed) grid points per period, 
respectively. These are solutions to equations (3.11) at each time step. Again, the 
dotted curve on this figure is the 200 time step per period result of figure (4.3b). We 
note that we get results identical to this control case, using as few as 2 time steps per 
period. 

We present two final figures associated with the linear simulations. Figure (4.6) 
shows Real(a?) versus time for the simulation with period fitting using 2 time steps 
per period, and demonstrates the effect of the mean flow evolution on the period 
of oscillation. In figure (4.7) we give a direct comparison between the growth rate 
predicted by OS theory at each time (dashed), versus the computed growth rate from 
the linear simulation (solid), again using the 2 time steps per period case. We see that 
the evolution has a significant destabilizing effect. 
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m(cj) 



0.28 0.32 0.36 0.40 


time 

Figure 4.6 Linear Sinnulation with period fitting. Duration: 20 
periods. 2 time steps per period. Period of 
oscillation as a function of time. 



Figure 4.7 Comparison of Orr— Sommerfeld growth rate computed 
locally at each time step (dashed) versus growth 
rate computed from linear simulation which allows 
growth of the boundary layer (solid). 
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4.3 Nonlinear Simulations 

For all the nonlinear simulations, we use the initial condition generated by the 2 nd 
order correction for the j = 1 mode, and its conjugate for the j = - 1 mode. However, 
we adjust the amplitude to test for the presence of nonlinear effects. All higher order 
modes are initially taken to have 0 values, as is the mean flow correction u$. 

For the first four sets of simulation results, we used 17 modes ( j = —8, -7,..., 8) 
and we fixed the peak amplitude of the initial profile for the fundamental mode, given 
in figure (4.2), to be 0.000001 (figure 4.8), 0.002 (figure 4.9), 0.005 (figure 4.10) and 
0.008 (figure 4.11). We employed the period fitting algorithm and used 40 time steps 
per one period of the fundamental mode. At each time step, the set of nonlinear 
equations given by (3.18) and (3.19), with boundary conditions (3.16) and (3.17) are 
solved. Each of these figures contains two graphs: the (a) part gives the growth rate 
of the fundamental mode, and in the (b) part, we present the energy curves for modes 
j = 1 to j = 8 where we define energy as ln(| v /F^|), and Ej is defined for the j th mode 
by equation (2.34). The strong oscillations exhibited m the higher order modes in 
figures (4.10) and (4.11) result from strongly nonlinear behaviour in the initial stages 
of the simulation, since all higher order modes are initially set to 0. This effect is 
amplified by aliasing errors. 

On each of these figures, the dotted curve represents the results from the linear 
simulation which uses the same initial conditions. For the 0.000001 case, the linear 
results are obtained over the entire time range (dotted curve sits on the solution for 
the fundamental mode). However, even when the amplitude is as small as 0.002, the 
nonlinear effects can be seen within just a few periods. 

To test the effects of even higher order terms, we repeat the 0.008 case above, using 
33 modes (j = —16,— 15,. ..,16). Again we use 40 time steps per one period of the 
fundamental mode. These results are presented in figure (4.12). In the (b) part of this 
figuie, we again show only the modes j = 1 through j = 8. A reduction in oscillations 
occurs due to the elimination of aliasing errors; however, the startup effects can still 
be seen. 
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Filially, to verify that the period-fitting algorithm and nonlinear algebraic solver are 
working successfully for these simulations, we include a simulation where explicit 2 nd 
order Adams-Bashforth temporal discretization is used to treat the nonlinear terms, 
and period fitting is not used. At each time step, we solve the linear set of equations 
(3.12) through (3.17). In this case, 400 time steps per one period of the fundamental 
period are used, and we include 16 modes. Figure (4.13) compares the growth rate of 
this simulation (20 periods in duration) with the first 20 periods of figure (4.11). It is 
evident that these algorithms are performing effectively. 
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Figure 4.8 Noniineor Simulation with period fitting. Duration: 50 
periods. 17 modes. Initial amplitude of fundamental 
mode: 0 000001. (a) Growth rate of the fundamental 

mode. (b) Energy curves. Solid curves in (b) correspond 
to the first and second harmonics; dashed curve is the 
mean flow correction. Dotted curve in both figures is 
the linear result. 





4 9 Nonlinear Simulation with period fitting. Duration. 50 
periods. 17 modes. Initial amplitude of fundamental 
mode: 0.002. (a) Growth rate of the fundamental mode 

(b) Energy curves. Dashed curve in (b) is the mean 
flow correction. Dotted curve in both figures is the 
linear result. 


E 





Figure 4.10 Nonlinear Simulation with period fitting. Duration: 50 
periods. 17 modes. Initial amplitude of fundamental 
mode: 0.005. (a) Growth rate of the fundamental mode 

(b) Energy curves. Dashed curve in (b) is the mean 
flow correction. Dotted curve in both figures is the 
linear result. 





Fiaure 4 11 Nonlinear Simulation with period fitting. Duration. 50 
periods 17 modes. Initial amplitude of fundamental 
mode: 0.008. (a) Growth rate of the fundamental mode 

(b) Energy curves. Dashed curve in (b) is the mean 
flow correction. Dotted curve in both figures is the 
linear result. 





Figure 4.12 Nonlinear Simulation with period fitting. Duration: 50 
periods. 03 modes. Initial amplitude of fundamental 
mode: 0.008. (a) Growth rate of the fundamental mode 

(b) Energy curves. Dashed curve in (b) is the mean 
flow correction. 





Figure 4. 1 3 Nonlinear Simulation without period fitting. Explicit 
treatment of nonlinear terms. Duration: 20 periods. 

17 modes. Initial amplitude of fundamental mode: 0.008 
Growth rate of the fundamental mode is shown. Dotted 
curve is simulation associated with figure (4.11) - 
obscured. 




§5 Conclusions 

In this article we have endeavoured to describe how Tollmien-Schlichting waves 
may develop in a temporally evolving flow. We have shown, as is to be expected, that 
above a certain critical amplitude their presence will lead the situation to undergo 
some kind of transition. We have not tried to describe the flow within this regime, we 

have purely tried to demonstrate a possible mechanism which may lead to this state 
of affairs. 

The main theme of our article is to develojD reliable numerical techniques which 
allow us to solve the nonlinear partial differential equations that can arise in this kind 
of problem. The methods that have been used have been drawn from the evolution of 
so-called parabolising numerical techniques, both for instability work (see Bertolotti 
(1991)) and for full Navier-Stokes simulations (El-Hady (1989)). In our work we have 
the distinct advantage that our underlying physical problem evolves with time rather 
than a downstream coordinate, which renders the governing equations parabolic rather 
than elliptic as in the spatially evolving boundary layer case. This means that it is 
not necessary to change the characteristics of the equations to effect a cheap numerical 
solution, rather the methods can be exploited to provide more accurate solutions whilst 
using far less resolution, for instance by period fitting. In the PSE work, the spatially 
evolving wavenumber is taken to vary linearly with the downstream coordinate, so that 
the governing equations are no longer elliptic. If the next order terms were retained, 
the equations would revert to their previous state, whereas in our problem we are free 
to retain as high order variations as we choose. The assumption we choose to make is 
that we can fit a periodicity and the remainder of the evolution is on a slow scale. 

In deriving initial conditions for our calculations we start with the Orr-Sommerfeld 
results, obtained by Otto (1994). As mentioned earlier we are aware that these results 
are fatally flawed, since they do not allow for the evolution of the undisturbed flow, 
but nevertheless they provide an initial guess for the modes’ form and frequency. We 
then consider an eigenproblem in which we not only solve for the mode and complex 
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frequency, but also the temporal derivative of both quantities. In the example con- 
sidered, this produces a fairly large change (rendering a mode that Orr-Sommerfeld 
predicts as stable to be unstable). This can be extended so that we ensure the second 
derivatives of the quantities are correct, producing minor modifications to the modes 
characteristics and eliminating just about all the visible transients from the growth 
rate curves (Figure 4.3). Unfortunately one problem that we have not resolved is the 
transients due to the start of the nonlinear modes; at the present we start our modes 
with zero amplitude (for the non-fundamental). We believe that this is what may be 
leading to the oscillations in the energy curves in figures (4.10b), (4.11b) and (4.12b). 

Our simulations for the linear problem seem to produce the expected results; that 
is, the flow can support Tollmien-Schlichting waves for a finite time interval. In the 
simulation associated with figure (4.7), it appears that the evolution of the mean flow 
has a significant destabilising effect, but it should be noted that the OS analysis, to 
which we are making comparisons, is suspect when applied to this evolutionary flow 
situation: in Section 4.1, we point out that the OS analysis predicts a stable mode 
where the corrected cases suggest instability (we appreciate that this mode is close to 
the neutral boundary). It should also be noted that even changing the quantity which 
represents the energy of the mode can significantly alter the results, as seen in Smith 
(1979a). Perhaps the main message to be taken from this part of the study is that one 
should consider the corrected Orr-Sommerfeld equation rather than the parallel one 
as the extra work is more than rewarded by the enhanced accuracy. 

The nonlinear simulations are as would be expected, in that above a certain initial 
amplitude of the fundamental, after a period of growth, the modes will induce a 
transition. The development of this set of methods has enabled us to study these 
modes, and we are currently attempting to take this calculation further, perhaps on 
into the turbulent regime. It is likely that our modal calculations will have moie 
success than the spatial calculations since our equations are actually truly parabolic. 

We hope to extend this work into the compressible domain, so that a study of the 
unstable characteristics in the case of subsonic, transonic and supersonic plate speeds 
may be performed, and the noise radiation at high speeds predicted. 
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